#include <vector>
#include "CUDA_full_weighting.h"
#include "CUDA_vec.h"
#include "CUDA_interpolation.h"
#include "CUDA_JacobiRelax.h"
#include "CUDA_Poisson.h"
#include "ConjugateGradient.h"
#include "utility.h"
#include <stdio.h>

#include "MultigridSolver.h"

#define DEBUG

using namespace std;

typedef CUDA_Full_Weighting RestrictionOp;
typedef CUDA_Interpolation InterpolationOp;
typedef CUDA_JacobiRelax RelaxOp;
typedef CUDA_Poisson Matrix;
typedef CUDA_vec<float> Vector;

int
main(int argc, char * argv[])
{
	int p = 10;
	int N = (1 << p) - 1;

	float P = 12;
	Matrix A(true); //use shared memory
	//Level 0 data
	Vector x1(N * N); Vector b1(N * N);
	x1.set_constant(P);
	print_CUDA_err_msg("Before multiply");
	A(x1, b1);
	print_CUDA_err_msg("Before zeroing");
	x1.set_constant(0);

	print_CUDA_err_msg("Before norming");
	printf("xnorm: %f bnorm: %f\n", norm(x1), norm(b1));

	print_CUDA_err_msg("After norming");
	MultigridSolver mg;
	mg(A, x1, b1, 1e-1);

	Vector true_x(N*N);
	true_x.set_constant(P);
	true_x.AXPBY(1,-1,x1);
	printf("Norm diff: %f\n", norm(true_x));

	return 0;
}
